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<^ Abstract 

Directional estimation is a common problem in many tracking applications. Traditional niters 
such as the Kalman filter perform poorly because they fail to take the periodic nature of the 
problem into account. We present a recursive filter for directional data based on the Bingham 
distribution in two dimensions. The proposed filter can be applied to circular filtering problems 
with 180 degree symmetry, i.e., rotations by 180 degrees cannot be distinguished. It is easily 
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implemented using standard numerical techniques and suitable for real-time applications. The 
presented approach is extensible to quaternions, which allow tracking arbitrary three-dimensional 
orientations. We evaluate our filter in a challenging scenario and compare it to a traditional Kalman 
£> " filtering approach. 
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1. Introduction 
"xl" \ 

Many estimation problems involve the task of estimating angular values. These problems include, 
but are not limited to, estimating the pose or orientation of objects. For example, tracking cars, 
ships, or airplanes may involve estimation of their current orientation or heading. Furthermore, 
many applications in the area of robotics or augmented reality depend on reliable estimation of 
the pose of certain objects. When estimating the orientation of two-way roads or relative angles of 
two unlabeled targets, the estimation task reduces to estimating an axis. This can be thought of 
as estimation of a directionless orientation or estimation with 180° symmetry. All these estimation 
problems share the need for processing angular or directional data, which differs in many ways from 
the classical Euclidean setting. First, periodicity needs to be taken into account. This is especially 
important for measurement updates around 0, respectively 2n. Second, directional quantities do 
not lie in a vector space. Thus, there is no equivalent to a classical linear model, as there are no 
linear mappings. 
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In many current applications, even simple estimation problems involving angular data are of- 
ten considered as traditional linear or nonlinear estimation problems and handled with classical 
techniques such as the Kalman Filter [l|], the extended Kalman Filter (EKF), or the unscented 
Kalman Filter (UKF) [2J|. In a circular setting, most traditional approaches to filtering suffer from 
assuming a Gaussian probability density at a certain point. They fail to take into account the 
periodic nature of the problem and assume a linear vector space instead of a curved manifold. This 
shortcoming can cause poor results, in particular when the angular uncertainty is large. In certain 
cases, the filter may even diverge. 

Classical strategies to avoid these problems in an angular setting involve an "intelligent" reposi- 
tioning of measurements or even discarding certain undesired measurements. Sometimes, nonlinear 
equality constraints have to be fulfilled, for example unit length of a vector, which makes it neces- 
sary to inflate the covariance [3J. There are also approaches, that use operators on a manifold to 
provide a local approximation of a vector space [J|. While these approaches yield feasible results, 
they still suffer from ignoring the true geometry of circular data within their probabilistic mod- 
els, which are usually based on assuming a normally distributed noise. This assumption is often 
motivated by the Central Limit Theorem, i.e., the limit distribution of a normalized sum of i.i.d. 
random variables with finite variance is normally distributed [5(. For angular data, this is not the 
case. Choosing a circular distribution for describing uncertainty offers possibly better results. 

In this paper, we consider the use of the Bingham distribution [6| for recursive estimation of 
orientation. The Bingham distribution is defined on the hypersphere of arbitrary dimension, so it 
can be applied to problems of different dimensionality. Here, we focus on the two-dimensional case 
and apply our results to axis estimation. To the best of our knowledge, this is the first published 
attempt to create a recursive filter based on the Bingham distribution. 

The presented methods can also be applied to the four-dimensional case, which would allow the 
representation of unit quaternions. Unit quaternions could then be used to estimate the full 3D 
orientation of an object. It is well known that Quaternions avoid the singularities present in other 
representations such as Euler angles. Their only downsides are the fact that they must remain 
normalized and the property that there are two quaternions for every rotation (q and — q). Both of 
these issues can elegantly be overcome by use of the Bingham distribution, since it is by definition 
restricted to the hypersphere and is 180° symmetric. 

This paper is structured as follows. First, we present an overview of previous work in the area 
of directional statistics and angular estimation (Sec. [2]). Then, we introduce our key idea in Sec. El 
In Sec. HI we give a detailed introduction to the Bingham distribution and we derive the necessary 
operations, which we will need to create a recursive Bingham filter. Based on these prerequisites, we 
introduce our filter in Sec. We have carried out an evaluation in simulations, which is presented 
in Sec. El Finally, we conclude this work in Sec. \7\ 

2. Related Work 

Directional statistics is a subdiscipline of statistics, which focuses on dealing with directional data. 
Classical results in directional statistics are summed up in the books by Mardia and Jupp [7( and 
by Jammalamadaka and Sengupta [a]. Directional statistics differs from traditional statistics by 
the fact that random variables located on manifolds (for example the circle or the sphere) are 
considered rather than random variables located in vector spaces (typically K ). 



There is a broad range of research for investigating the 2D orientation, e. g., the work by Krindis 
et al. [9(. A recursive filter based on the von Mises distribution for estimating the orientation on 



distributions was presented in 



the SO{2) was presented in [10 1. L ater, a nonlinear filter based on von Mises and wrapped normal 
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In 1974, Bingham proposed his distribution in [6(. Further work on the Bingham distribution 



has been done by Kent [12| as well as Jupp and Mardia 131 ] . So f ar, there have only been a few 
applications of the Bingham distribution, for example in geology [14|. In 2011, Glover used the 
Bingham distribution for a Monte Carlo based pose estimation [15|. Glover also released a library 



called libbingham Iff] that includes implementations of some of the methods discussed in Sec. |U 



It should be noted that our implementation is not based on libbingham. 

3. Key Idea of the Bingham Filter 

The goal of this paper is the derivation of a recursive filter based on the Bingham distribution. 
Rather than relying on the traditional Gaussian distribution, we chose to represent all occurring 
probability densities as Bingham. The Bingham distribution is defined on the hypersphere and 
is antipodally symmetric, which makes it interesting for applications in angular estimation with 
inherent 180° symmetry and for problems where 180° symmetry occurs as a result of parameteri- 
zation, e.g., in the case of quaternions. Although we restrict ourselves to the two-dimensional case 
in this paper, we would like to emphasize that most of presented methods are easily generalized to 
higher dimensions. 

In order to derive a recursive filter, we need to be able to perform two operations. First, we 
need to calculate the predicted state at the next time step from the current state and the system 
noise affecting the state. In a traditional estimation problem in M. d with additive noise, this involves 
a convolution with the noise density. We provide a suitable analogue on the hypersphere, which 
we call composition. Since Bingham distributions are not closed under compositions, we present 
an approximate solution to this problem, which is based on matching covariance matrices. 

Second, we need to perform a Bayes update. As usual, this requires the multiplication of the 
prior density with the likelihood density. We prove that Bingham distributions are closed under 
multiplication and show how to obtain the posterior density. 

4. Bingham Distribution 

The Bingham distribution appears naturally when a d-dimensional normal random vector x with 



E(x) = is conditioned on ||x|| = 1 [17|j. In the following, we will introduce the Bingham dis- 
tribution and derive the formulas for multiplication of two Bingham probability density functions. 
Furthermore, we will present a method for computing the composition of two Bingham-distributed 
random variables, which is analogous to the addition of real random variables. 

4-1. Probability Density Function 

Definition 1. Let S^-i = {x € K rf : ||x|| = 1} C K d be the unit hypersphere in IR d . The probability 
density function (pdf) 
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Figure 1: Bingham pdf with M = 12x2 and Z = diag(— 8,0) as a 3D plot. This corresponds to a 
standard deviation of 16°. 



of a Bingham distribution \d^ is given by 
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exp(x T M Z M T x) , 



where M € M. is an 

Z = diag(zi, . . . z d -i,0) € R dxd with z 1 <_ 

normalization constant. 



orthogonal matrix (NL~M. T = M^M = Idxd) describing the orientation, 
< Zd-i < is the concentration matrix, and F is a 



As Bingham showed, adding a multiple of the identity matrix Idxd to Z does not change the 
distribution. Thus, we conveniently force the last entry of Z to be zero. Because it is possible to 
swap columns of M and the according diagonal entries in Z without changing the distribution, we 
can enforce z± < ■ ■ ■ < z^—i- This representation allows us to obtain the mode of the distribution 
very easily by taking the last column of M. 

The pdf is antipodally symmetric, i.e., f(x) = /(— x) holds for all x 6 S^—i- Consequently, the 
Bingham distribution is invariant to rotations by 180°. Examples of the pdf for two dimensions 
{d = 2) are shown in Fig. Q] and Fig. [5J The Bingham distribution is very similar to a Gaussian if 
and only if the uncertainty is small. This can be seen in Fig. [3l which shows the Kullback-Leibler 
divergence between a Bingham pdf and a corresponding Gaussian pdf. 

4-2. Normalization Constant 

The normalization constant can be calculated with the help of the hyper geometric function of 
a matrix argument 18|, [l9|, 12CJ] . It is given by 



\Sd-i\ ■ iF\ ( 2'2' Z 




Figure 2: Bingham pdf with M = I2 X 2 for different values of Z = diag(zi,0) and x - 
(cos(#),sin((9)) . These values for z\ correspond to standard deviations of approximately 36 c 
16°, and 6° respectively. 



where |<Sd_i| is the surface area of the d-sphere and ii*i(-, •, •) is the hyper geometric function of 
matrix argument. In the two-dimensional case [d = 2), this reduces to 
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so it_is sufficient to compute the hypergeometric function of a scalar argument, which is described 
in 



21] 



4-3. Multiplication 

For two given Bingham densities, we want to obtain their product. This product is used for 
Bayesian inference involving Bingham distributions. The result presented below yields a convenient 
way to calculate the product of Bingham distributions. 

Lemma 1. Bingham distributions are closed under multiplication with renormalization. 

Proof. Consider two Bingham distributions 

fi{x) = Fi • exp(x T Mi Zi Mjx) 

and 

f 2 {x) = F 2 ■ exp(x T M 2 Z 2 M^x) . 
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Figure 3: Kullback-Leibler divergence on the interval [0, n] between a Bingham pdf with M = 
I2x2 ; Z = diag(zi,0) and a Gaussian pdf with equal mode and standard deviation. For small 
uncertainties (z\ < — 15, which corresponds to a standard deviation of about 11°), the Gaussian 
and Bingham distributions are almost indistinguishable. However, for large uncertainties, the 
Gaussian approximation becomes quite poor. 



Then 
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fi{x) ■ / 2 (x) = F X F 2 ■ exp(x J (MiZiMf + M 2 Z 2 M£ ) x 



oc F ■ exp(x T M Z M T x) 

with F as the new normalization constant after renormalization, M are the unit eigenvectors of 
C, D has the eigenvalues of C on the diagonal (sorted in ascending order) and Z = D — D^I^xd 
where D^d refers to the bottom right entry of D, i. e., the largest eigenvalue. □ 

4-4- Estimation of Parameters 

Estimating parameters for the Bingham distribution is not only motivated by the need to 
estimate noise parameters from samples. It also plays a crucial role in the prediction process when 
computing the composition of two Bingham random vectors. This procedure is based on matching 
covariance matrix. Be aware that although the Bingham distribution is only defined on S^—i, we 
can still compute its covariance in Mr. Thus, we will present both the computation of the covariance 
matrix of a Bingham distributed random vector and the computation of parameters for a Bingham 
distribution with a given covariance (which could originate from an arbitrary distribution on the 
hypersphere) . 

The maximum likelihood estimate for the parameters (M, Z) of a Bingham distribution can be 
obtained as described in [g]. M can be obtained as the matrix of eigenvectors of the covariance 



S with eigenvalues uj\ < C02- In other words, M can be found as the eigendecomposition of 
S = M • diag(wi,u;2) ■ M T . To calculate Z, the equations 
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l-Fi(£,Ml) 

have to be solved under the constraint 22 = 0, which is justified by the argumentation above and 
used to simplify the computation. This operation is performed numerically. 

Conversely, for a given a Bingham distribution (M, Z), the covariance matrix can be calculated 
according to 

S = M • diag(wi, w 2 ) • M T 

__ ,. (\ dF 1 8F\ A;rT 
= M • diag — -— , — — - • M 1 

1 dF 



£^~-M(:,;)M(:,i 






where M(:,i) refers to the z-tli column of M (161 ] . Thus, any Bingham distribution is uniquely 
defined by its covariance matrix and vice versa. The following Lemma simplifies the computation 
of partial derivatives of a confluent hyper geometric function of a 2 x 2 matrix argument, which is 
used in computation of the covariance matrix as derived above. 

Lemma 2. For d = 2, the partial derivatives 

can be reduced to hypergeometric functions of scalar argument. 



Proof. See Appendix A □ 



4-5. Composition 

Now, we want to derive the composition of Bingham distributed random variables, which is 
the directional analogue to adding random variables. This operation can, for example, be used 
to disturb an uncertain Bingham-distributed system state with Bingham-distributed noise, similar 
to using a convolution to disturb a probability distribution on R with additive noise. First, we 
define a composition of individual points on the hypersphere S^-i, which we then use to derive the 
composition of Bingham distributions. 

The composition of two Bingham distributions depends on the interpretation of the unit vectors, 
for example as complex numbers or quaternions. We assume that a composition function 

© : S d _i x S d _i -> S d _i 

is given. The function © has to be compatible with 180° degree symmetry, i.e., 

±(x © y) = ±((-x) © y) = ±(x © (-1O) = ±((-s) © {-v)) 



for all x, y G Sd-i- Furthermore, we require the quotient (Sd_i/{±1}, ©) to have an algebraic group 
structure. This guarantees associativity, the existence of an identity element, and the existence of 
inverse elements. 

In the complex case, we interpret S\ C M 2 as unit vectors in C, where the first dimension is 
the real part and the second dimension the imaginary part. In this interpretation, the Bingham 
distributions can be understood as a distribution on a subset of the complex plane, namely the 
unit circle. 

Definition 2. The composition function © is defined to be complex multiplication, i.e., 

xi\ fyi\ = fxiyi - x 2 y 2 
x 2 ) \V2j \x1y2 + x 2 yi 

analogous to 

(xi + ix 2 ) ■ (yi + iy 2 ) = (xiyi - x 2 y 2 ) + i(xiV2 + x 2 yi) . 

Since we only consider unit vectors, the composition © is equivalent to adding the angles of 
both complex numbers when they are represented in polar form. The identity element is ±1 and 
the inverse element for (xi,x 2 ) is the complex conjugate ±(xi, —x 2 ). 

Unfortunately, the Bingham distribution is not closed under this kind of composition. That 
is, the resulting random vector is not Bingham distributed. Thus, we propose a technique to 
approximate a Bingham distribution to the composed random vector. The composition of two 
Bingham distributions /a and /b is calculated by considering the composition of their covariance 
matrices A, B and estimating the parameters of /c based on the resulting covariance matrix. 
Composition of covariance matrices can be derived from the composition of random vectors. 

Lemma 3. Let /a and /b be Bingham distributions with covariance matrices 

an ai 2 \ fbn b\ 2 

and rs = 

* 022/ V * 022 

respectively. Let x,y £ Si C I 2 be independent random vectors distributed according to /a and /b- 
Then the covariance 

'en cia j . = Cov (^ ey ) 

. * C22/ 
of the composition is given by 

en =011611 - 2ai2&i2 + 022622 , 

C12 =011612 - ai 2 6 2 2 + 012611 - 022612 

C22 =011622 + 2012612 + 022611 . 



Proof. See Appendix B □ 



Based on C, maximum likelihood estimation is used to obtain the parameters M and Z of the 
uniquely defined Bingham distribution with covariance C as described above. This computation 
can be done in an efficient way, because the solution of the equation involving the hypergeometric 
function is the only part which is not given in closed form. This does not present a limitation to 
the proposed algorithm, because there are many efficient ways for the computation of the confluent 
hypergeometric function of a scalar argument 22J, |23 |. 



Algorithm 1: Prediction 



Input: estimate M|,Z|, noise M^,Z^ 

Output: prediction M^ +1 , Z^ +1 

/* calculate covariance matrices A, B */ 



a\ 



A = Eti^gM|(:,i)M|(:, 

B = Eti^f£Mn:,z)M»(:,i) T ; 

/* calculate C according to Lemma [3] */ 

en = anbn - 2ai 2 h 2 + a 2 2&22! 

C\1 = au&12 - 022^12 - 012&22 + «12&li; 
C 2 2 = ail^22 + 2ai 2 6l2 + 022&11; 

'en Ci 2 N 

v Cl2 C 2 2, 

/* calculate M| +1 , Z|! +1 based on C */ 

Mg +1 ,ZP +1 <-MLE(C); 



5. Filter Implementation 

The techniques presented in the preceding section can be applied to derive a filter based on the 
Bingham distribution. The system model is given by 

x k+1 = s* © m , 

where w k is Bingham distributed noise. The measurement model is given by 

z k = x k ev k , 

where v k is Bingham distributed noise and x k is an uncertain Bingham distributed system state. 
Intuitively, this means that both system and measurement model are the identity disturbed by 
Bingham distributed noise. Note that w k and v k can include a constant offset. For example 
w k could include a known angular velocity. Alternatively, to avoid dealing with biased noise 
distributions, a rotation may be applied to x k first and unbiased noise added subsequently. 

The predicted and estimated distributions at time k are described by their parameter matrices 
(M|!, 7i P k ) and (M|, Z|) respectively. The noise distributions at time k are described by (M]f, Zjj?) 
and(M£,Z£). 

5.1. Prediction 

The prediction can be calculated according to 

(M p k+1 , Z p k+1 ) = composition((M|„ Z|), (M£, Z£)) , 

which uses the previously introduced composition operation to disturb the estimate with the system 
noise. 



Algorithm 2: Update 



Input: prediction M£, Z£, noise M^,Z^, measurement z_ k 



Output: estimate M|,Z 



e 
ki k 



/* rotate noise according to measurement */ 

M^(IeM^); 

/* multiply with prior distribution */ 

(M%, Z%) <- multiply((M, Z£)), (M*, Z*)); 



5.2. E/pdate 

Given a measurement 2, we can calculate the updated distribution according to Bayes' rule 

/(M fc , ZfclD = c • /(l|M fe , Z fc ) • /(M fc , Z fc ) 

with some normalization constant c, which yields the update procedure 

(M£, Z|) = multiply((M, Z e fc ), (M£, Z|)) 

with M = (z MJ), where a indicates the complex conjugate of a and is evaluated for each 
column of MY. 

6. Evaluation 

The proposed filter was evaluated in simulations. In this section, all angles are given in radians 
unless specified differently. 

For comparison, we implemented a one-dimensional Kalman filter [l[. A traditional one- 
dimensional Kalman filter has two issues when confronted with our situation. First, it does not 
take the circular nature of the problem into account. Second, it does not handle 180° symmetry. 
We can circumvent both issues by restricting the estimate x k according to < x^ < ir and by 
shifting the measurement, so that \xk — Zk\ < § is satisfied. 

In our example, we consider the estimation of an axis in robotics. This could be the axis of a 
symmetric rotor blade or any other robotic joint with 180° symmetry. We use the initial estimate 
with mode (0, 1) T 

mi=(\ °V z* ( ~ x ° 



10 " \0 l) 1 
the system noise with mode (1,0) 

n?=r? !V z- 



" l 1 0/ ' k \ 



and the measurement noise with mode (1,0) 
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MJ= i • Z t 



M TV- ~ 3 ° 

i o > Zfc -\ o o 
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(a) Ground truth and estimate. 




40 60 

time step 



(b) Angular error. 
Figure 4: An example run of Bingham and Kalman filter. 

The true initial state is given by (1,0) T , i.e., the initial estimate with mode (0, 1) T is very poor. 
The initial estimate for the Kalman filter is given by 

Xq = atan2(mode(Mo)) = atan2(l,0) 



7T 



and the noise means are 

H™ = n v k = a tan2(0, 1) = . 

The covariance matrices for the Kalman filter are obtained by sampling the Bingham noise 
parameters and calculating the empirical covariance from the samples. This yields 

Cg = 0.5956, C% = 0.0027, C v k = 0.2836 , 

which is equivalent to standard deviations of 44° for the first time step, 3° for the system noise and 
30° for the measurement noise. 

We simulate the system for a duration of fe max = 100 time steps. An example run is depicted 
in Fig. [H In addition to the mode of the estimate, we plot the 95% confidence interval, which is 
equivalent to the 2a bounds in the case of the Kalman filter. 
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For evaluation, we consider the angular RMSE which is given by 



\ 



k 



'"max 

-Efoo 2 

max k=l 



with angular error 

e k = min(Z(xf lc ,mode(M^)),7r - 4a£ ue ,mode(M|)) 

at time step k. Obviously, < e k < § holds, which is consistent with our assumption of 180° 
symmetry. 

The presented results are based on 1000 Monte Carlo runs. Even though our filter is computa- 
tionally more demanding than a Kalman filter, it is still fast enough for real-time applications. On 
a standard laptop, our non-optimized implementation in MATLAB needs approximately 60 ms for 
one time step (prediction and update), which could be significantly improved by a faster evaluation 
of the hypergeometric function. In Fig. [6l we plot the error of our filter against the error of the 
Kalman filter for all runs. The proposed filter outperforms the Kalman filter in most cases, which 
is also true for the mean angular error in every time step as shown in Fig. [5j In particular, the 
significantly faster rate of convergence of the proposed filter is evident. This superiority is due 
to the reasons listed in the introduction. The use of the Gaussian distribution, which does not 
consider the problem geometry leads to suboptimal results compared to the proposed approach 
based on the Bingham distribution. 

In Fig. O we also show a comparison with a filter based on the wrapped normal distribution 



(denoted WN), which we previously published in [ll[ and modified for the 180° case. The angular 



error of both filters is almost indistinguishable. However, unlike the proposed filter based on 
the Bingham distribution, the previously published filter cannot easily be generalized to higher 
dimensions. 

7. Conclusion 

We have presented a recursive filter based on the Bingham distribution. It can be applied to 
circular estimation problems with 180° symmetry. Our simulations have shown the superiority of 
the presented approach compared to the traditional solution of modifying a Kalman filter for the 
circular setting. 

Future work will focus on recursive 3D pose estimation using Bingham distribution. This can be 
achieved by applying the presented methods in the four-dimensional case for estimating quaternions. 
Open challenges include an efficient estimator of the Bingham parameters based on available data. 
This makes an efficient evaluation of the confluent hypergeometric function necessary. 
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40 60 

time step 



(a) Time steps 1 < k < 100. 




(b) Time steps 20 < k < 100. 



Figure 5: Results of 1000 Monte Carlo runs. We show the mean error at every time step across all 
runs for the proposed Bingham filter, a Kalman [l( filter and a filter based on the wrapped normal 
(WN) distribution ll[. Because the initial error is large as a result of the poor initial estimate, we 
show two plots of different time intervals. 



Appendix A. Proof of Lemma [2J 

We use the identities 

1 ■, (zi ° 



and 



This yields 



l-ri 



F, -,1, 



z 2 



exp(z 2 ) • i-Fi ( 2' 1 ' 21 ~ Z2 



— 1 F 1 {a,b,z) = ^ 1 F l (a + l,b+l,z) . 



g F = A F (1 1 fa ° 

0*i |5 d _i| a^i 1 x V 2 ' ' v° ^ 

d (\ 

= exp(z 2 )^— i-Fi I -, 1, zl - zl 

= exp(z 2 )-iF 1 I -, 2, zl - z2 
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Figure 6: Results of 1000 Monte Carlo runs. Each sample represents one run. Samples below the 
red line indicate that the proposed filter has performed better, samples above the red line indicate 
that the Kalman filter has performed better. 



and 



g I = A_ F (I i (* o 

dz 2 \S d _ x \ 3z 2 l l \2' '\0 z 2 

= ~dz ( ex P( Z2 ) li?1 ( 2 ,1,Zl ~ Z2 

= exp(z 2 )iF 1 ( 2>l>zi ~ z 2 \ +exp(z 2 )— 1 F 1 l-,l,zi - z 2 

= exp(z 2 )iFi (-,l,z 1 ~ z 2j -exp(z 2 )-iFi ( t^' 2 ' 21 - z 2 
= exp(z 2 )l l-Pi ( -,l,zi -z 2 j - g i-^i ( 2 ,2,Zl ~ Z2 



D 



14 



Appendix B. Proof of Lemma [3], 

The covariance of the composition 
C = Cov(x y) 

=coJ( xm - x * y A) 

\\xiy 2 + x 2 y 1 J J 

= /Var(a;iyi - x 2 y 2 ) Cov(x 1 y 1 - x 2 y 2 , x x y 2 + x 2 yi) 
V * Yai(xiy 2 + x 2 yi) 

can be obtained by calculating the matrix entries individually. For the first entry we get 

Cll =Var(xiyi - x 2 y 2 ) 

= E((xiyi - x 2 y 2 f) - {E(xiyi - x 2 y 2 )f 

=E(x\yl - 2x iyi x 2 y 2 + x\y\) - (E(sij/i) - E(x 2 y 2 )) 2 (B.l) 

= E{x\) E{y\) - 2 E( Xl x 2 ) E(y m ) + E(^) E(y 2 2 ) (B.2) 

-(E(x 1 )E(y 1 )-E(x 2 )E(y 2 )) 2 (B.3) 



=011611 - 2cii 2 bi 2 + a 2 2&22- 

We use independence of x and y in (jB.ip , linearity of the expectation value in ()B.2[) and symmetry 
of the Bingham in (|B.3[) . Analogously we calculate 

C22 =011^22 + 2a\ 2 b\ 2 + a22&n ■ 
The off-diagonal entry can be calculated similarly 

C12 = CovOiyi - x 2 y 2 ,x 1 y 2 + x 2 y\) 

= E((xiy\ - x 2 y 2 ) ■ (xiy 2 + x 2 yi)) - E(x i y 1 - x 2 y 2 ) ■ E{xiy 2 + x 2 yi) 
= E(xly ± y 2 - x x x 2 y\ + x x x 2 y\ - x\y x y 2 ) - (E(cci) E(yi) - E(x 2 ) E(y 2 )) 

■ (E( Xl )E(y 2 ) + E(x 2 )E( yi )) 
=anb 12 - a\ 2 b 22 + ai 2 6n - 022^12 • 

D 
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